\chapter{Setting up, running and restarting simulations}

\section{System generation}\label{sec:init}

\subsection{Pure water systems}\label{sec:initWat}

The simplest simulation in \brahms involves ``pure'' water systems, that is, systems containing only water molecules. The relevant parameters in {\tt brahms.md} are discussed in the following.

\paragraph{{\tt initUcell}} Initial coordinates are set on a face-centered cubic (FCC) lattice following a procedure taken from Rapaport~\cite{rapa}[p.67]. 
The lattice contains 4 sites per unit cell, and the number of unit cells is set through {\tt initUcell}, which is a 3-element integer parameter.
For example, by setting {\tt initUcell 6 6 6} \brahms generates a system containing $4\times(6\times6\times6)=864$ sites arranged on the FCC lattice.

\paragraph{\tt nWaters} The parameter {\tt nWaters} should be set to the total number of water sites, hence to $4*initUcell.x*initUcell.y*initUcell.z$.

\paragraph{\tt nSites} The parameter {\tt nSites} should be set to the total number of sites, which in this case corresponds to {\tt nWaters}.

\subsection{Single-component lipid bilayers}\label{sec:initLip}

To generate a bilayer system from scratch, \brahms requires a few parameters to be set in the input file {\tt brahms.md}.

\paragraph{\tt nLipids}
First of all, the number of lipids {\tt nLipids} must be chosen. The algorithm that generates the initial configuration requires {\tt nLipids} to be set such that the quantity $\sqrt{(nLipids/2)}$ is an integer. For example, by choosing 512 {\tt nLipids}, the quantity $\sqrt{(nLipids/2)}$ is 
equal to 16; in this case, \brahms will create two ``mirrored'' monolayers each constructed with 16*16 lipids arranged on a regular grid.
It follows that acceptable values for {\tt nLipids} are for instance: 32, 50, 72, 98, 128, 162, 200, 242, 288, 338, 392, 450, 512, \dots, 1058, \dots, 5000, etc.

\paragraph{\tt nWaters}
The number of water molecules hydrating the bilayer is determined through the parameter {\tt nWaters}. This should ideally be based on experimental data. 
In particular, experimental studies often describe the hydration level of the bilayer in terms of the parameter ``n$_W$'', that is, the number of water molecules per lipid.
In this case, the target value for {\tt nWaters} should be $n_W\times nLipids$. In some cases, experimental studies report the hydration level $H$, defined as $H = 100 \times weight_{water} / weight_{total}$.
As for {\tt nLipids}, also {\tt nWaters} cannot be set to any value, because of the underlying algorithm that adds water to the system. However, there is enough flexibility
to obtain a value which is going to be ``close enough'' to the desired one - see next paragraph.

\paragraph{{\tt initHalfCellWat}}
In particular, the generation of water is controlled by the parameter {\tt initHalfCellWat}, which ``initializes'' water in each of the two ``halves'' of the system,
above and below the bilayer. This parameter actually contains 3 elements, which are 3 integer numbers controlling the generation of the initial water ``lattice''.
For example, setting {\tt initHalfCellWat} to 10 10 5 creates two water slabs, above and below the bilayer, each characterized by 5 ``planes'' of water molecules, 
with each plane containing $10\times 10$ water molecules in the corresponding xy plane; therefore, in this example, the total number of water molecules will 
be $2\times (10\times 10 \times5)$, where the factor~2 accounts for the generation of 2 slabs (above and below the bilayer). By varying the three integers in 
{\tt initHalfCellWat} it is possible to generate a number of water molecules to approximate the desired value. Example values for the input parameter {\tt initHalfCellWat} are 
collected in table~\ref{tab:genWat}.
\begin{table}
\begin{center}
%\vspace{14pt}
\begin{tabular}{c|c|c} % put @{} if you want to eliminate horizontal space between columns
{\tt nWaters} & {\tt initHalfCellWat}& notes\\ 
\hline
1014 &  13 13 3& suitable for 32 DOPCs\\
1024 & 16 16 2 & suitable for 128 DOPEs\\
1536 & 16 16 3& suitable for 128 gel-DSPCs\\
2304 & 24 24 2 & suitable for 288 DOPEs\\
2400 & 20 20 3 & suitable for 200 gel-DSPCs\\
2500 & 25 25 2 &\\
3072 & 16 16 6 &\\
3240 & 18 18 5 & \\
3610 & 19 19 5 &\\
3872 & 22 22 4 &\\
4000 & 20 20 5 &\\
4046 & 17 17 7 & \\
4056 & 26 26 3 &\\
4096 & 32 32 2 & also 16 16 8\\
4232 & 23 23 4 & suitable for 128 DOPCs\\
10890& 33 33 5 &\\
12544& 56 56 2 & suitable for 1058 gel-DSPCs\\
16384 & 32 32 8 &\\
16428 & 37 37 6 &\\
16810 & 41 41 2 & suitable for 512 DOPCs\\
16854 & 53 53 3 &\\
16928 & 46 46 4 &\\
33800 & 65 65 4 &\\
34848 & 66 66 4 & suitable for 1058 DOPCs\\
\end{tabular}
\caption[Examples for {\tt initHalfCellWat}]{Examples for {\tt initHalfCellWat} (nWaters=2*initHalfCellWat.x*initHalfCellWat.y*initHalfCellWat.z).}
\label{tab:genWat}  
\end{center}
\end{table}
{\em NOTE:} as you can see from the table, the third value of {\tt initHalfCellWat} is normally set to $2-8$, irrespectively of the total amount of water. This is because the thickness of the water phase ``above and below'' the membrane is roughly constant (for standard hydration levels) no matter how large is the bilayer surface. 

\paragraph{\tt nSites} The parameter {\tt nSites} should be set to the total number of sites, hence $nSites = nLipids * nSitesPerLipid + nWaters$, 
where $nSitesPerLipid$ is 15 in the case of the current ELBA models for DOPC, DSPC and DOPE. 

\paragraph{\tt region} This input parameter is a 3-element vector or real numbers which sets the lengths of the edges of the simulation box along the x, y and z directions.
The first two values (x and y) represent the edges of the bilayer plane, by convention. These can be set considering the experimental area per lipid $A_L$:
$$ region.x = region.y =\sqrt{A_L * nLipids / 2} $$
where the division by 2 accounts for the fact that each monolayer is made of $nLipids/2$ lipids.
The last element, {\tt region.z}, sets the length of the simulation box along the z direction, by convention perpendicular to the bilayer plane. The total volume 
of the simulation region (the ``box'') is:
$$ region.x \times region.y \times region.z = V_L \times nLipids + V_W \times nWaters $$
with $V_L$ the experimental volume per lipid and $V_W$ the experimental volume per water molecule (0.03\,nm$^3$ at room temperature). Hence we can write:
$$ region.z = \frac{V_L \times nLipids + V_W \times nWaters}{region.x \times region.y} $$


 \section{Restarting simulations}\label{sec:restarting}

Restarting a simulation is an important procedure in Brahms. In fact, most simulations are typically not started from scratch, but restarted from some previously equilibrated configuration. \brahms implements the checkpoint/restart method described in Rapaport~\cite{rapa}(pp.~500-504).

To enable checkpointing, the {\tt doCheckpoint} input parameters must be set to 1; also, {\tt stepCheckpoint} must be given a value representing the checkpointing frequency (as number of steps). \brahms will then periodically write the checkpoint files {\tt mdXXchecka.data}, {\tt mdXXcheckb.data} and {\tt mdXXcklast.data}, where {\tt XX} corresponds to the input parameter {\tt runId}. 

When checkpointing is enabled, at the beginning of the run \brahms automatically looks in the directory for the checkpoint files. If they are present, the corresponding simulation is restarted. If there are no checkpoint files, a new run begins and the related checkpoint files are generated.

When restarting a simulation, by default \brahms retains the step count from the corresponding previous run. This means that the restarted simulation will run for ({\tt stepLimit} - previousStepCount) steps.\footnote{So for example if {\tt stepLimit} $<$ previousStepCount, \brahms will terminate immediately after restarting.} Sometimes it can be convenient to reset the step count; this can be easily done by setting {\tt resetTime} to 1. In this case, the restarted simulation will run for {\tt stepLimit} steps.

Note: the checkpoint files {\tt mdXXchecka.data}, {\tt mdXXcheckb.data} and {\tt mdXXcklast.data} are binary files, so they might not be portable when architectures are (very) different. 

\section{Self-assembly simulations}

This section describes how to set up self-assembly simulations. If an initial ``random'' system is not available, it can be generated by ``disassembling'' a bilayer at high temperature.

\subsection{``Disassembling'' run}
Protocol to ``scramble'' a bilayer into a random mixture for subsequent self-assembly simulation:
\begin{enumerate}
\item  reduce timestep to 10 fs
\item  increase temperature to 1000-2000\,\textcelsius
\item  turn thermostat on, switch barostat off (NVT conditions)
\item  set cutoff radii to 0.9\,nm (optional)
\item  set 'hBondXX' parameters to 1
\item  set charges and dipoles to 0 
\item  run for as long as it takes to obtain a random configuration; regularly visualize the system to check progress
\end{enumerate}


\subsection{Equilibration run}

\begin{enumerate}
\item  restore desired temperature while keeping the barostat switched off
\item  restore desired force field parameters and cutoff radii
\item  run with small timestep (e.g., 1\,fs) to pre-equilibrate (remember to set {\tt tauT} accordingly - see~\textsection\ref{sec:brahmsmd}).
\item  restore normal timestep (e.g., 15\,fs) and run for as long as required
\item  switch on the barostat ({\tt applyBarostat 1}) under isotropic control ({\tt flexBox 0}) and run for as long as required
\end{enumerate}

\subsection{Self-assembly run}

Starting from a random configuration, a self-assembly run can be carried out like a normal job. Two options can be considered to control the pressure:
\begin{itemize}
\item isotropic (more robust): {\tt flexBox 0}
\item anisotropic (maybe more realistic, but potentially leading to instabilities): {\tt flexBox 1} and {\tt keepSquare 0}
\end{itemize}
Snapshots from a typical self-assembly simulation are reported in Figure~\ref{fig:self}.
\begin{figure}[ht]
\centering
\includegraphics[scale=0.3]{runningSimulations/self.eps}
\caption[Bilayer self-assembly]{Bilayer self-assembly. Snapshots from a self-assembly simulation of a system containing 128 DOPC lipids and 4900 water molecules. Simulation length was $\approx50$\,ns.}
\label{fig:self}
\end{figure}
